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ABSTRACT 


Extension to the Dynamic Modeling of the 
Large Angle Magnetic Suspension Test Fixture 
Lucas E. Foster 
Old Dominion University 


The Large Angle Magnetic Suspension Test Fixture (LAMSTF) is a 
laboratory scale proof-of-concept system. The configuration is unique in that the 
electromagnets are mounted in a circular planar array. A mathematical model of 
the system had previously been developed, but was shown to have inaccuracies. 
These inaccuracies showed up in the step responses. Eddy currents were found 
to be the major cause of the modeling errors. In the original system, eddy 
currents existed in the aluminum baseplate, iron cores, and the sensor support 
frame. An attempt to include the eddy current dynamics in the system model is 
presented. The dynamics of a dummy sensor ring were added to the system. 
Adding the eddy current dynamics to the simulation improves the way it 
compares to the actual experiment. Also presented is a new method of 
determining the yaw angle of the suspended element. From the coil currents the 
yaw angle can be determined and the controller can be updated to suspend at 
the new current. This method has been used to demonstrate a 360 degree yaw 
angle rotation. 
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CHAPTER ONE 


Introduction to Magnetic Suspension Technology 

Magnetic suspension is a relatively new science. While it is currently 
being used in some industrial and scientific applications, its use will increase 
dramatically in the future when the science is perfected. 

Magnetic suspension applications can be divided into two groups 
according to the dimension of the air-gap, which is the distance between the 
object being suspended and the suspending magnets, as compared to the overall 
scale of the system. These groups are large-gap and small-gap. The difference in 
gap size is important because the assumptions made on a given problem will be 
different. An example of an assumption that can be used in small-gap system 
that is invalid for a large-gap system is that the fields are uniform across the air 
gap and have no gradients [1]. This assumption makes small-gap systems easier 
to analyze. 

In a typical small-gap suspension system the gap can be as small as a 
fraction of a millimeter. Magnetic bearings and very accurate small object 
positioning systems are examples of small-gap systems. Magnetic bearings serve 
the same purpose as regular bearings, which is to hold a rotating shaft in place, 
yet they have several advantages over conventional bearings. The parts in the 
bearing do not touch so wear is minimal and friction is eliminated. These 
characteristics translate into longer life with less down-time, and higher 
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efficiency. Magnetic bearings are currently being used in a wide range of 
applications such as in pump stations along pipelines[2]. They are of great 
advantage because they may be in a remote location that makes performing 
routine maintenance costly. Magnetic bearings should not require any 
maintenance, and they have no reason to ever wear out. 

One known application of small gap magnetic levitation is the 
magnetically levitated train [3,4]. They have several advantages over typical 
trains that run on rails. They are faster, quieter, more efficient, and smoother. 
These characteristics are all the result of removing the friction associated with 
wheels and the rumblings of all the mechanical parts in an ordinary train. They 
are typically faster because of two reasons. The first reason is that the train only 
has to overcome the force of air resistance. The second reason is in the design of 
the track. The track can be designed in such a way as to be safer with respect to 
derailing as compared to conventional railroad tracks, and the maglev track cam 
be sloped to make banked turns allowing higher speeds than would be wise in a 
conventional train. Another advantage of a magnetically levitated train is that 
the train is powered magnetically allowing the power to come from an energized 
track. This cuts the weight of the train by eliminating the engine, a substantial 
amount of weight, making it even more efficient. 

Other types of small gap suspension systems include accurate moving 
systems that are used in the manufacture of integrated circuits [5]. In these 
applications being able to move a circuit board a few Angstroms with accuracy is 
important. Another advantage in using magnetic levitation is the ability to 
isolate the board from machine vibration 

Large-gap magnetic suspension systems are extremely varied in their 
applications. For a system to be classified as a large gap system, the gap may 
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be as little as 0.1 meters, for a small system, or much larger. Some of the uses 
for large gap systems include pointing, vibration isolation, and wind tunnel 
model support applications. 

Pointing systems have a suspended element that can be controlled about 
different axes to aim a payload [6]. Pointing can be used by telescopes mounted 
on space stations that axe in orbit. The telescopes need to be able to track a star 
across the sky, and by suspending the telescope it can track the star smoothly 
and accurately. 

Vibration isolation systems are used in microgravity experiments where 
vibrations from equipment not related to the experiment would invalidate the 
experiment [7]. Vibration can be transmitted to the experiment through the 
support structure and umbilicals, or the experiment itself may create vibration. 
Using levitation the experiment no longer touches the vibrating platform and the 
other sources of vibrations can be actively eliminated. 

Magnetic Suspension and Balance Systems (MSBS) are used in wind 
tunnels. In most ordinary wind tunnels, the models are suspended by fixing 
them to the walls of the tunnel through the use of a model mount called a sting. 
The problem is that the sting can alter the airflow around the model making the 
experiment invalid [8]. The requirement of the model to be altered to accept the 
sting may also make the experiment invalid. However, by suspending the model 
magnetically, the need for a sting is eliminated along with a major source of 
errors. The amount of force on the model can be derived from measuring the 
currents going through the coils and subtracting them from the currents used 
when the model is being suspended without any external force. A computer, 
which is most likely controlling the wind tunnel anyway, can determine how 
much force, such as drag or lift, is on the model. 
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There are some magnetic bearings and large gap suspension systems that 
axe configured so they axe stable without any type of active controls [9]. 

However, most practical applications of magnetic suspension technology require 
an active controller. An active controller measures the position of the object 
being suspended and develops an error signal indicating how fax from the datum 
position it is. Based on this error signal the controller changes the currents in 
the coils to develop forces on the suspended object to push or pull it back to its 
datum position. Depending on the configuration the controller may have to 
update the currents fifty times per second or more. 

Magnetic suspension can be achieved in a number of different ways. The 
method used by magnetic bearings relies on the fact that ferromagnetic material 
is attracted to a magnet. The shaft is made from a ferromagnetic material like 
iron or steel, and it is surrounded by electromagnets and position sensors at the 
location of the bearing. The suspension is achieved by attractive forces pulling 
the shaft back to the center. 

Other suspension systems, the Large Angle Magnetic Suspension Test 
Fixture (LAMSTF) in particular, rely on magnetic interaction between a 
permanent magnet that is the suspended element and the electromagnetic coils. 
LAMSTF is a set of five electromagnets mounted in a flat planar array and will 
be described in more detail in chapter two. The suspended element is a magnet 
with its north and south poles at the ends, and the suspension is maintained by 
producing an electromagnetic south pole under the south pole and an 
electromagnetic north pole under the north pole of the suspended element. The 
coils are used to change the fields and the gradients of the field in subtle ways to 
keep the suspended element from falling out of suspension. For LAMSTF the 
forces are all repulsive, but this is not true in general. In most other applications 
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electromagnets axe placed overhead to provide attractive forces as well as 
repulsive forces. 

Chapter three develops the analytical model of the system, but this model 
has been shown to have inaccuracies causing the dynamics of the suspended 
element to be different from what is expected. A large portion of the system 
errors are due to the eddy currents in the system. Chapter four attempts to 
correct the errors by developing a way to mathematically describe the eddy 
currents and add them to the system. 
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CHAPTER TWO 


The Large Angle Magnetic Suspension Test Fixture 

The LAMSTF was constructed as a proof-of-concept system for a unique 
hardware configuration. The hardware configuration is unique for several 
reasons. The coil arrangement, five coils of the same size set in a horizontal 
circular arrangement, suspends the element using only repulsive forces. This 
form of suspension is called levitation, but the terms suspension and levitation 
will be used interchangeably. Most other large gap systems use a combination of 
attractive and repulsive forces to suspend. The LAMSTF system configuration is 
also unique from a controls standpoint because there is a lack of natural 
decoupling. In most other systems the coils are arranged so that each coil affects 
only one or two degrees-of-freedom, but in LAMSTF each coil has an effect on all 
of the degrees-of-freedom. This aspect makes control more difficult because the 
controller must decouple the different degrees-of-freedom. More information 
covering all aspects of LAMSTF including the controls, modeling and potential 
variations of the configuration can be found in [10]-[14]. 

The hardware used for the LAMSTF can be broken down into three 
different sections. The first section is a description of the plant. The plant 
includes the magnets, amplifiers, and the suspended element. The plant consists 
of the elements of the system that determine the dynamics. The other two 
systems are the position feedback sensors and the computer that is used for 
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control. These systems are interconnected by means of voltages. The sensors 
produce six voltages that are related to the position of the suspended element. 
These electrical voltages are used by the computer to find a set of five voltages 
that are related to the desired currents in the five electromagnets. 

2.1 The Plant Description 

The system consists of five electromagnets arranged in a planar array in 
the form of a circle under the suspended element. The coils are arranged 
symmetrically at a radius of 13.77cm. This configuration is shown in Figure 2.1. 
The coils have an outer diameter of 15 cm, and an inner diameter of 8.5 cm. 

They are wet wound using epoxy on bakelite forms with 509 turns of AW G 10 
wire. In order to keep them stable while they are running a coating of epoxy has 
also been applied to the outside of the coils. In the center of each coil is an iron 
core. The purpose of the core is to boost the output of the coils so they can 
suspend the element using reasonable currents. The iron core is as long as the 
electromagnet with an outer diameter of 7.5 cm. The core is not laminated, 
which is common practice to reduce eddy currents, and this poses problems that 
will be explained in the section on eddy current modeling. 

The electromagnets are driven by amplifiers that take a voltage between 
± 10 volts. They amplify the voltage by -3, and they have a output current 
between ± 30 amperes. For safety reasons the amplifiers are connected to 
temperature probes located inside the coils. If a coil reaches 160 degrees all of 
the amplifiers axe shut off automatically. This is to reduce the risk of damaging 
the coils, and it takes about 30 minutes of solid operation to produce that 
temperature. 
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Under the current configuration, the system can control the model in five 
degrees- of -freedom, and they are x, y, z, translations and pitch and yaw 
rotations as shown in Figure 2.2. The system has no control in the roll degree-of- 
freedom. As a result of this, during suspension the element will rock freely until 
it is stopped by physically constraining the rolling motion. Another system is 
under development that has more coils for the purpose of controlling the 6*^ 
degree-of-freedom. This system is also a prototype proof-of-concept system that 
uses repulsive forces. 

2.2 The Position Sensors 

The sensors that axe currently used on LAMSTF axe light blocking position 
sensors. They have two components, a light emitting component, and a 
receiving component. The emitter discharges an infrared laser that has a 
rect an gular cross-section area that is about .5x.5 inches. This beam has a 
uniform amount of light spread out across its cross section. The detector 
translates the amount of light reaching it to a voltage. When the suspended 
element is floating at its center or zero position it will block half of the light 
beam. Then, when the model moves around, the amount of light will change 
and hence the amount of voltage coming from the detector will change. This 
voltage typically changes from one for a fully blocked beam to five volts for a 
clear beam. 

In an area of high magnetic fields that are constantly changing, such as 
when the element is being suspended, the noise induced in a long voltage run can 
become significant. This problem is solved by placing a current loop in the 
voltage line. This setup allows a large portion of the run to be a current based 
signal which is less prone to noise caused by magnetic field interference. Where 
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a voltage signal is required, the wires are twisted together to reduce the noise. 

In LAMSTF, measurements from six sensors are used in locating the 
suspended element along the five degrees-of-freedom. The sensors are arranged 
in a circle in a planar array above the magnets (see Figure 2.3). Sensors 1 and 2 
are used to find z position and pitch. Sensors 3, 4, 5, and 6 are coupled to find 
the x, y, and yaw position. After the voltages are brought into the computer via 
the analog-to-digital converter board, they are multiplied by the sensor 
decomposition matrix S2P in the following equation. 


Position = S2P * Sensor Voltages 


S2P= 


0 

0 

.353 

-.353 

-.353 

.353 

0 

0 

-.353 

.353 

-.353 

.353 

.5 

.5 

0 

0 

0 

0 

1 

1 

0 

0 

0 

0 

0 

0 

-.5 

-.5 

.5 

.5 


After the multiplication is performed, there are five numbers which 
represent the five degrees-of-freedom x, y, z, pitch, and yaw. It is these five 
numbers which are used by the controller as inputs. 

Currently the sensors are the limiting factor in large angle control of the 
model. The sensor range depends on the size of the model. The sensors measure 
one diameter of the model in translation. This is because if the model is high 
enough to allow light to show under the model, as in Figure 2.4, the sensors will 
indicate the wrong position, and the controller will not respond correctly, 
dropping the suspended element out of suspension. 
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2.3 The Computer’s Role as a Controller 


The computer’s responsibility is to read in the sensor voltages and 
produce current output voltages. As the model moves about, the current in the 
coils must be constantly updated in order to maintain stable suspension. The 
rate at which the computer updates the currents can vary, but for a typical run 
500 Hz is used. One of the additional goals of the project is to develop a 
controller that can run at 40 Hz, but this is difficult because it is very slow 
compared to the fastest unstable mode of oscillation. The computer uses three 
boards: a counter-timer board, an analog-to-digital converter board, and a 
digital- to- analog converter board. 

For the controller routine to work it has to be run at a constant and 
known frequency. However, the program used to control the element has a 
graphical interface to the user which is very slow, and the amount of time it 
takes to run varies because it updates different amounts of the screen depending 
on the different settings. In order to make sure that the control routine is being 
run at a constant rate the program uses a hardware interrupt. A hardware 
interrupt is used by other things like the printer and mouse. In all, there are 256 
different interrupts possible in the computer. Only a small portion of them, 
eight to be exact, can be used by a programmer. The purpose of the counter- 
timer board is to generate the interrupt, and the board can be set up from the 
program in an initialization routine. The program sets the board to generate a 
square wave signal at the desired frequency. This signal is fed into the interrupt 
pin on the board. Whenever a rising edge is sensed by the interrupt controller 
chip it stops whatever the computer is doing and runs the interrupt service 
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routine. For the purposes of the magnetic suspension the interrupt service 
routine is set to be the controller routine. After the controller has been run a 
signal is sent to the interrupt controller chip to tell it that the routine is finished 
and it can return to whatever the computer was doing before the interrupt 
occurred. This method allows the graphics, which are slow, to proceed at their 
own rate while the controller routine interrupts it at the known frequency needed 
by the controller. 

The other boards are used to translate the voltages used by the amplifiers 
and the voltages from the sensors to numbers that can be used by the computer. 
The analog to digital converter takes the 0 to 10 volts output by the sensors and 
converts it to a number from 0 to 4095. This number is an integer that is 
proportional to the voltage. The digital-to-analog board repeats the process 
starting with a number and converting it to a voltage that is used by the 
amplifiers. 

The interrupt service routine is the controller. It uses the six sensor 
voltages to produce the five positions. The controller dynamically compensates 
the five position signals with lead or lag to generate the five coil currents. For 
the controller used in this thesis, the coil currents are updated 500 times a 
second. 
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CHAPTER THREE 


General Dynamic Modeling 

An accurate mathematical model of the system is necessary to create a 
control law that will be successful in suspending the element. Utilizing the 
mathematical model, a potential controller response is simulated to access the 
closed-loop stability and performance. An iterative design process is used to 
eliminate any deficiencies uncovered in the simulation. The system that is being 
modeled is the plant as described in Chapter 2. 

The general procedure for developing the system model starts with the 
basic equations. These basic equations are linearized and put into a useful form 
to be used in simulation. 

3.1 The Coordinates 


There are two different coordinate systems used to make the analysis 
simpler. Both of the coordinate systems are typical right handed orthogonal 
systems. The first system is called the inertial coordinates. They are fixed with 
respect to the coils, and the origin is at the zero position of the suspended 
element. The second system is called core coordinates, and they are allowed to 
move around as the suspended element moves. The x axis is along the length of 
the suspended element with the y axis horizontal and the z axis vertical. 
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When the model does not show any translation or rotation from its zero 
position the two coordinate systems are as shown in Figure 3.1. Then as the 
suspended element moves around the core coordinates track the motion. The 
core coordinates are used to find the torques and forces on the model that result 
from the different fields and gradients. The inertial coordinates are used to find 
the fields and gradients. In actual use the field and gradients will be expressed 
in inertial coordinates and transformed to core coordinates. Transposing from 
one coordinate system to the other is a simple matter of developing a 
transformation matrix, based on the location of the model. In general operation, 
the coordinate system will be very close together due to the fact that all of the 
step responses will be kept to small angles and translations. In this situation, 
transformation will not be necessary. However, one degree of freedom will 
undergo large angular displacement and it is the yaw angle. The system has 
demonstrated a yaw angle of 360 degrees. If the z axis remains coincident with 
the inertial z axis then the transformation matrix will look like: 


T m = 


cos(0 z ) 

-sin(0 z ) 

0 


sin(0 z ) 

cos(0 z ) 

0 


0 

0 

1 


where 0 Z is the yaw angle. This transformation matrix will transform field 
vectors and velocity vectors. Transforming a derivative between the inertial 
coordinate system and the core coordinate system is not as easy. This is because 
the core coordinates are allowed to accelerate. The equation for transforming the 
force and torques, the derivatives of linear and angular momentum, is as follows: 

F = m c + fi x m c V 
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where V is the velocity of the core , m c is the mass of the core, and ft is the 
angular rotation rate of the moving core coordinates with respect to the inertial 
coordinates and it is expressed in matrix form as 


[ft] 


0 -ft- ft 
z 

ft- 0 -ft. 



0 


This is the typical procedure for developing equations in an accelerating 
coordinate system such as the core coordinates. For the LAMSTF system the 
accelerations will, under most circumstances, be small enough to be neglected. 
This assumption is made to greatly simplify the development of the governing 
equation. 

3.2 The Basic Equations 


The torques and forces due to the magnet fields are: 

T = J v (M x B) dV 
F = J v (M x V) B dV 


These two equations are used to solve for the linear accelerations and the 
angular acceleration. T and F are the torque and force vectors, M is the 
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magnetic vector, and B is the field vector. The integrals are evaluated over the 

volume of the core. The first assumption is that the equations can be rewritten 
as: 

T = V (M x B) 

F = V (M • V) B 


This assumption treats the field like a constant and is possible because 
throughout the core the fields are symmetric functions about the suspension 
point. The modeling errors introduced by the previous equations are constant 
values, and they are taken care of in the VM term, which is not capable of being 
measured with accuracy. 

At this point the fields axe in the inertial coordinates and the magnetic 
vector is in core coordinates. This is corrected by using the transformation 
matrix, as described earlier, in the equations as follows: 


T = V M [T m ] B* 

F = V [T m ] [3B] (T4 1 M 

These equations are used to solve for the linear and rotational 
accelerations as shown below. 


« = l£j VM [T m ] B 
^ i V |T m ] OB] [14 1 M 


17 



These equation will be used to develop the A and B matrices of. the linear state 
space representation of the system. The equation for velocity in x will be 
developed here. All of the other degrees of freedom are developed in the same 
fashion with a more complete development found in Reference [16]. The 
equation for V x is as follows. 


Vx — m c V (® xx + 2 0 Z Bxy - 2 Oy B xz j + Oy g 


A Taylor series expansion is performed on V x resulting in: 


V x = V x |„ + ^|„«X 


where X is the state variable in the state space representation. 


X = (0 y 0 Z X y z fty fl z V x V y V Z ) T 
A linear approximation for the change in velocity along x is then: 



SX 


When taking the partial of the equation, the field terms B are functions of 
x, y, z, and the coil currents. Two other assumptions are made to get to the 
next step. The first is that the Oy and 0 Z are small angular displacements. This 
assumption is reasonable due to the fact that the sensors can only sense a small 
angle. The second assumption is that products of small numbers, such as OySx, 
are negligible. The resulting 6V X is the equation actually modeled in the state 
space equation below. 
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6X = A *SX + B*6I 


Where 6l is the input currents that are used to create the fields. The B is 
not to be confused with the fields and gradients B which is denoted as a vector 
or with a direction subscript. The A matrix in the previous equation is: 


0 

0 

0 

0 

n 

Dxz 


0 

B X z 

Byz 

Bzz. 

0 

0 

0 

0 

n 

0 


n x 

Byz 

"Byy 

Byz 

0 

0 

0 

0 

0 


xx'j 

2D X y 

Bxxx 

Bxxv 

B X X7, 

0 

0 

0 

0 

0 

B\z 

(B 

yy'B xx ) 

B xyx 

Bxyy 

B X yz 

0 

0 

0 

0 

0 


By* 

B*7.X 

Bxzy 

Bx77 

1 

0 

0 

n 

0 

0 


n 

n 

0 

n 

0 

1 

0 

n 

0 

n 


n 

0 

n 

n 

0 

0 

1 

n 

0 

n 


n 

n 

0 

n 

0 

0 

0 

l 

n 

0 


0 

n 

0 

0 

0 

0 

n 

0 

l 

0 


0 

0 

0 

n 


Where W contains the constants from the equations along its diagonal. 
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The B matrix is made up of terms that are the inputs of the system. However, 
before the equation can be used the current from the coils must be introduced 
into the field inputs, B. This is done by finding the K values in the following 
equation. 

B = [K] f 


K is m units of gauss per amp. A K value exists for each different field 
Bx> By, and B z , and all of the gradients. This group of K values will become the 
B matrix in the state space equations. Following the same logic as before, a 
Taylor series expansion for the velocity can be written as: 


V x = V x + 


dl 


SI 
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Where SV X is: 


The SV X is the final equation that is in the B matrix. The B matrix is 

also multiplied by the W matrix that holds the constants. The B matrix is as 
follows: 

- R zi -K 22 -K Z3 -K Z4 

K y 2 *S K y 4 

Kxx i Kxx 2 k xx 3 k xx< 

Kx yi K *y 2 Kx y 3 K xy 4 

K xzj K XZ2 K XZj K XZi 

0 
0 
0 
0 
0 

The system that has been developed is a pure analytical model. This 
model has inaccuracies in it, and these inaccuracies show up in the step response 
of pitch. The damping of the step response in pitch is less in experiment than 
the simulated results predicted. A method of improvement for this problem is 
presented in Chapter 4. 




21 



3.3 System Analysis 


Now that the model has been developed the next step is to analyze the 
system to find out what modes exist and which modes are stable or unstable. 
The values for K are calculated from OPERA, which is a software package used 
to numerically evaluate magnetic fields. The method that OPERA uses to 
calculate the fields is described in Chapter 5. That chapter will also cover how 
the field gradients are calculated. The numeric values for the A and B matrices 
for the LAMSTF are as follows. 


A = W*\ 
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0 
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B = W* 


4.071el 

0 

2.3010e-l 

0 

-2.88G9e-l 


4.071cl 

9.527cl 

-1.5935e-l 

1.208e-l 

8.91Glc-2 


4.071el 

5.8899cl 

8.1846c-2 

-2.0525e-l 

2.3356c! 


4.071cl 
-5.8899cl 
8.1846c-2 
2.0525c- 1 
2.335Gc-l 


4.071cl 

-9.5278cl 

-1.5930-1 

-1.2GSe-l 

-8.91Glc-2 


0 

0 

0 

0 

0 


W = 


4.3308c5 0 

0 4.3308e5 •. 

1 1.0602e2 

1.0G02c2 

1.0G02c2 


1 



1 0 
0 1 


Using these matrices the eigenvalues for the plant, model are: 



Eigenvalues = 


59.2596+0i 

-59.2596+0i 

0+7.9717i 

0-7.9717i 

0+.9556i 

0-.9556i 

58.2942+0i 

-58.2942+0i 

9.7762+0i 

-9.7762+0i 


There axe five different pairs of eigenvalues corresponding to five different 
modes of motion. Two of the motions which lie on the ju> axis are neutrally 
stable. The other three modes are very unstable and they will throw the 
suspended element out of suspension. Table 3.1 has all of the eigenvalues and 
what degrees-of-freedom they correspond to. The two most unstable modes are 
called the compass needle modes which include yaw and x. The reason for the 
highly unstable behavior is because when the model is suspending, there is a 
large B x field component that is in the opposite sign as the magnetic vector in 
the suspended element. The B x component is not used to control any degree-of- 
freedom, but it exist when the other field and gradient are created. The desire of 
the suspended element is to align itself with the B x field causing the high 
frequency of the compass needle modes. Figure 3.2 has a diagram of the 

different modes of the suspended element as they relate to the motions of the 
model. 
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3.4 Initial conditions 


The actual currents through the coils will be the offset current and the 
control current. Finding the initial condition is important for knowing what the 
offset current through the coils should be. The initial currents can be found 
using statics and the original two equations that the model was built upon. The 
only force on the model will be the force required to offset gravity. By setting 
all of the other forces and torques to zero the fields and gradients necessary to 
support the model can be found. 


By — B z — B xx — B X y — 0 


Bxz — 


_ me g 


V M-v 


Using the field equation form before: 


B = [K] I 


The initial currents will be: 


I'o = [K ] 1 B- 


The initial coil currents are listed in Table 3.2. The currents in coils 2 
and 5 are the same and the current in coil 3 and 4 are the same. This is due to 
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the symmetry of the circular planar configuration. Coils 1, 2, and 5 are all 
negative while 3 and 4 are positive. This is due to the fact that the magnetic 
vector of the suspended element is along the x axis. 
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Mode 

Eigenvalue 

Stability 

Degrees-of-freedom 

1 

± 59.26 rad/s 

Unstable 

x, 9y (Axial, Pitch) 

2 

± 7.972 i rad/s 

Stable Oscillatory 

x, 9y (Axial, Pitch) 

3 

± 58.29 rad/s 

Unstable 

0 Z (Yaw) 

4 

± .0956 i rad/s 

Stable Oscillatory 

z (Vertical) 

5 

±9.776 rad/s 

Unstable 

y (Lateral) 


Table 3.1 The Eigenvalues of the System 


Coil 


Current (Amps) 


1 

2 

3 

4 

5 


-14.17 

-4.38 

11.47 

11.47 

-4.38 


Table 3.2 Coil Currents 
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Figure 3.1 The Axis System 



Vertical translation (stable) 


Lateral translation (unstable) 


Figure 3.2 Modes of Motion 
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CHAPTER FOUR 


Eddy Current Analysis 

The discrepancy in the step response in pitch was believed to be due to 
eddy currents. In the original design of the system, aluminum was used as a 
base plate and as the support structure for the sensors. Iron cores that are not 
laminated were used as the core of the electromagnets in order to produce 
stronger fields with less current. All of this metal permits eddy currents to flow 
which affects the system plant giving the model errors. 

Eddy currents are produced anytime a conducting material is placed in a 
fluctuating magnetic field. As the magnetic field changes, currents are developed 
in the metal that resist the change. The fields calculated for a given set of 
currents through the coils is accurate for a steady-state condition. However, 
when the current going through the coils is changing the field response is not 
immediate. The transient field at the suspension point is reduced by the eddy 
currents in the system. In the original system dynamic model the eddy currents 
were ignored. However, the model was sufficiently close to the actual system 
dynamics to allow an implemented controller design to suspend the model. 
Discrepancies in the transient behavior between the analytical model and the 
physical system still exist though, and correcting these deficiencies is the main 
focus of the research. By developing the eddy current analysis, the theoretical 
simulation will move closer to the actual responses. 
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In LAMSTF the eddy currents arise from 3 different places. It was 
suspected that the greatest effect came from the sensor support ring. This effect 
shows up most prominently when a step response in pitch is performed (see 
Figure 4.1). The aluminum base plate that the coils are mounted on has an eddy 
current associated with it, but it is believed that its effect is reasonably small. 

The iron cores are the most challenging to model accurately. The reason for this 
is explained later in the theoretical analysis. A figure of all the eddy currents is 
shown in Figure 4.2. 

In most other magnetic suspension systems, eddy currents are minimized 
by design. Eddy currents need a closed path to travel or they cannot exist. The 
effects of the eddy currents can be reduced by electrically cutting the eddy 
current loops. Loops can be cut by using laminated cores and non-conducting 
inserts in the frame. Eddy currents can be avoided entirely by replacing the 
metal with a nonconductive material such as wood or fiberglass. 

As was stated before, it was suspected that the most prominent eddy 
currents come from the sensor support frame. The sensors are mounted on a 
large aluminum ring that is placed over all of the magnets. The reason for pitch 
to be particularly affected is when a step response in pitch is performed, the 
amount of field passing through the ring has a large net change. Step responses 
in the other degrees-of-freedom are not as sensitive to the sensor ring eddy 
currents. This is due to the fact that when the step responses are performed the 
field is changed in such a way as to leave the net amount of field passing through 
the ring the same. For these reasons the step response in pitch will be used as a 
comparison between the original model and the updated model. 
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4.1 Theoretical Outline 


The eddy currents will be modeled as a transformer with the secondary 
coil shorted out. This method will give a good first-order approximation of most 
eddy currents as long as they satisfy a few assumptions. The first assumption is 
that the skin depth, the depth that the eddy currents penetrate into the 
material, is much greater than the thickness of the metal. This also implies 
that the eddy currents are uniformly distributed throughout the metal. The skin 
depth varies according to the frequency that the field is being driven. As the 
frequency is increased the skin depth decreases. The frequency that the skin 
depth is calculated at in order to check this assumption is the highest frequency 
that will affect the system; for the case of LAMSTF the highest unstable open 
loop pole is about 20 Hz, to allow for a safety factor 100 Hz is used. The 
equation for skin depth is as follows[15]: 


skin depth — y fj. 0 fi r c u 

where /x r is the permeability of the material, fi 0 is the permeability of free space 
which is a constant, a is the resistance of the material, and u> is the field 
frequency. 

The large skin depth assumption is used throughout the following 
theoretical calculations, but considering the case of the iron cores as the 
frequency is increased the assumption loses its validity. The model s 
inaccuracies at high frequency may cause the theory to be off slightly at the 
frequencies of interest. The reason for the deviation at high frequencies is that 
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as the frequency is increased, the magnetic fields and eddy currents are forced to 
the outside of the core leaving the middle of the core free of magnetic fields and 
currents. This makes the uniform current assumption invalid. As the frequency 
gets very high the currents are reduced to surface currents along the outside of 
the core. This reduction in the amount of material that the eddy currents are 
flowing through increases the resistance of the loop, one of the parameters that 
was assumed to be constant in the equation for the eddy current dynamics that 
will be developed. 

Figure 4.3 is used as a general reference for the analysis. All of the eddy 
currents use the same equation with different parameters as will be demonstrated 
later. 

For the driven coil, the voltage across the te rminals is: 


V = IR + L & + L : 




ii 

i dt 


+ L 


m 


fe2 

2 dt 


+ 


The first two terms deal with the primary coil itself without any eddy 
currents. Each of the other individual terms deal with one eddy current. L m is 
the mutual inductance between the n^ eddy current loop and the driving coil, 
and I 6n is the amount of current in the loop. The mutual inductance L mn is a 
constant. For each eddy current loop, the following equation can be used to 
eliminate I en from the previous equation. The following equation does not 
include the effect of the other eddy currents on the eddy current being 
considered. If it were desired the equation could be amended to reflect the other 
eddy current loops, but for the sake of simplicity the effects of eddy currents on 
eddy currents has been assumed to be negligible. 
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For the eddy current loop the following equation applies: 


0 - Ie, Re, + L e , “3T + L m $ 

where all terms noted with an e are the values in the eddy current loop itself. 
The last term is the mutual inductance and the current in the driving coil. By 
taking the Laplace transform of both equations and solving: 




1 

\ 

(R + Ls) - 

’ 



Re, + ke, s ^ 


Now: 


I - 1 

v (R +Ls) 


is the normal transfer function for the terminal characteristics of the primary coil 
without any eddy currents. The other terms axe the effects of the eddy currents. 

To show that this is correct a special case is considered. The special case 
comes about in two different ways; the first is when the frequency is allowed to 
go to zero, and the second is the case where the resistance becomes very small. 

To make these points clear, two substitutions are made: 


L — qLp 
c i 


and 



33 



The result of which is: 



( 


(R + Ls) - 


V 


P 2 a(L ei sy 

(Rej + Re! S ) 


When R Ci is allowed to go to zero or as s becomes very large: 



1 

R + Ls (l-£ 2 


)) 


This indicates that the secondary loop effectively ” turns off” part of the 
primary coil if the secondary is reactive(non-dissipative). 

The field at any point from a coil can be calculated by: 



where Kj is the field coefficient and q is the number of coils to be added 
together. Eddy currents are added in the same manner. 





p 

n=l 


K 


e n, Ie n 


Taking the Laplace transform for one coil and eddy currents: 


-P K e L m s 


B- = K.I 1 - V D a ' 
J J 1 n=l K j + W) 
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Where j is the direction of the field desired and p is the number of eddy currents 
that are to be included in the equation. 

Now that the form of the equation has been developed the problem is to 
find values for the variables of a given eddy current loop. These variables can be 
calculated or variants of them can be derived from experimental method. The 
program used to calculate the variables is called VF / GFUN which is part of 
OPERA and it is distributed by Vector Fields Inc. VF/GFUN is a static field 
calculation program, but it is sufficient to find the variables that are needed in 
the above equations. 

Finding the resistance of a current loop is a relative simple task that 
involves finding the area of the loop, its length, and its resistivity which is a 
property of the material. These properties are used in the following equation. 


R = 



Length 

Area 


) 


Where R is the resistance, Length is the circumference of the loop, Area is the 
area of the cross-section of the loop. 

For the purposes of finding the other parameters, the eddy current loops 
are treated like a small coil of the same size as the assumed loop size. The K e 
value can be found from the equation 


B z = K z I 


where B z is found from OPERA when one coil or loop is turned on to one amp. 
For self-inductances, the eddy current loop is turned on to one amp and the field 
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going through the loop is integrated. The equation governing the voltage in the 
loop is: 


V — T dl _ 

L dt “ dt 


Therefore: 

oo 

Jb z dA 

t d<?i> o 

L “ dT“T- T — 

With the current at one volt, the integral becomes the inductance. For mutual 
inductances the driven coil will be the actual coil in the system. 

The accuracy of the results obtained by this method varies greatly 
depending on the geometry of the eddy current. If the eddy current is of a 
known geometry such as the sensor ring, then the results will be reasonably 
accurate to use in the model. Unfortunately, when the geometry of the eddy 
current is not of an exact or known size then the accuracy of the results are 
degraded according to the amount of simplification necessary. For the plate the 
diameter of the eddy current loop is not known, but it is assumed to be 
approximately the size of the inner diameter of the driving coil. No assumption 
is required concerning the height of the eddy current coil because the height is 
equal to the thickness of the plate. For the iron core the true eddy current 
geometry changes with different frequencies with the thickness and the height 
being the unknowns. For the core, the calculations serve as a gross magnitude 
approximation for selecting the values from the experimental results. 

There are magnetic field software tools that use numerical solution 
techniques, similar to OPERA, that will find where the eddy currents are flowing 
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in a piece of metal. ELEKTRA is one such program which is also marketed by 
Vector Fields, the same company producing OPERA. 

The preceding section explained how the parameters are evaluated from 
pure mathematical modeling. The next section will explain how to find the 
parameters from the measured transfer functions. A comparison of the two 
methods will be at the end of the next section. 

4.2 Experimental Procedure 

Experimental determination of the various transfer functions is necessary, 
not only to validate the theoretical calculations, but to improve the variables in 
the equations that will be added to the system. The transfer function can be 
determined from the experimental Bode plots. 

In order to get the experimental Bode plots, a Schlumberger transfer 
function analysis machine, model number 1250 was used. This analyzer drives a 
system with a sine wave, and the output of the system is fed back to the 
analyzer. In doing this, one point on the Bode plot is created. By stepping the 
sine wave through a given range of values, the frequency region of interest on the 
Bode plot can be generated. This data is stored in the machine as a file, which is 
then transferred to the computer via a GPIB connection. MATLAB is a 
mathematical analysis program that will be used to analyze the data. 

In order to get the desired Bode plot from the system that is being tested, 
several instruments are used. The amplifiers used for the system are the same 
ones used for all of the tests. The sine wave is fed directly into the amplifier 
driving one coil. The configuration of the metcil and the driving coil is changed 
in the test apparatus to isolate the effects of one eddy current. The field is 
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measured at a point with a F.W.Bell model number 9903 gaussmeter. The 
meter gives the field in one direction at one point using a Hall-effect sensor. The 
output of the Bell gaussmeter is an analog signal that is proportional to the 
amount of field at the probe. This signal is fed back to the transfer function 
analysis machine to form the input to output dynamic model. 

The desired Bode plot has an input of coil current and an output of field. 
In reality, it is not possible to get the exact transfer function. Every piece of 
equipment adds its own unique transfer function to the final Bode plot. To deal 
with this problem, a baseline measurement is taken that includes all of the 
different transfer functions. This baseline is subtracted from all other data 
taken, leaving behind only the desired Bode plot. The baseline configuration is 
an air- cored coil with the field probe lined up along the axis of the coil. The 
probe was placed at the suspension height for all tests including the baseline. A 
plot of the baseline is shown in the Figure 4.4. 

There are three different configurations that will be calculated and 
compared to experimental values. These three configurations are the iron core, 
the aluminum base plate, and the sensor ring. The configurations used for 
calculating and testing purposes are designed to isolate and identify the effects of 
one eddy current loop at a time. Two of the configurations, the aluminum base 
plate and the iron core, are exact representations of the corresponding eddy 
currents in the real system. However, the sensor frame is a very complex shape 
with many possible paths for eddy currents to travel. This will be modeled in 
the experiments by a horizontal aluminum ring. A picture of the system with 
the dummy sensor ring in place is shown in Figure 4.5. It is believed that this 
ring has the greatest effect on the system, and it is the easiest eddy current loop 
to calculate. 
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Using the experimental Bode plot and knowing the general form of the 
transfer function, it is possible to determine the parameters. This is done in 
MATLAB using a trial and error method. The general form of the transfer 
function that is used to get the parameters is as follows: 



The variables a and b axe varied to get the resonant frequency and the 
phase shift at the resonant frequency correct. The variables for the experiment 
and the theory are given in appendix B. The a variable will shift the phase 
difference while the b variable will change the resonant frequency. Finding the 
variables will not give the same parameters that are in the theoretical outline, 
but they will give the transfer function that will be used. The figures comparing 
the transfer functions will determine the validity of the approximations invoked 
throughout the theory. For the case of the core, the assumption of large skin 
depth becomes inaccurate as the frequency becomes large. At large frequency, 
the phase of the core starts to roll off whereas the theory returns to zero. To 
take care of this problem a new theory will have to be developed. 

The transfer functions given by the test data can be used to update the 
mathematical model of the system. In Figure 4.6 comparison of the test data, 
the theoretical modeling, and the experimental approximation are shown. 

4.3 Comparison of Experimental Results 


The results of the previous section are presented in Figure 4.5 through 
4.8. The figures contain two plots, the first one is of magnitude vs. frequency 
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and the second one is the phase plotted against the frequency. Figure 4.6 shows 
the response of the eddy currents in the dummy sensor ring. This plot contains 
the curves for the experiment, the original theory, and the estimated theory. 

This plot is presented to show how close the theory is to the actual experiment. 
For the other plots the experimental and the estimated theoretical curves axe 
shown without the original theory. This is because the differences between 
experiment and original theory axe due to inaccuracies in the estimation in 
parameters such as resistance, and the shape of the eddy current. Figure 4.7 and 
Figure 4.8 axe the effects of the eddy currents in the aluminum plate and the 
iron core. The estimated theory and the experiment agree reasonably well in the 
case of the aluminum plate and the dummy sensor ring. However, in the case of 
the iron core the agreement between estimated theory and experiment falls apart 
at the higher frequencies. The theory involved is not sufficient to explain these 
differences, but another theory is being developed that will better approximate 
the core response. That theory is still in its development stage and is beyond the 
scope of this report. Figure 4.9 is a response of all of the eddy currents in the 
system. Most of the error that exists in the plots is due to the discrepancy in the 
iron core modeling. This is given as an example of how several different eddy 
currents can be added together. 

4.4 Incorporation into Dynamic Model 

In the current system the aluminum plate has been replaced with a wood 
base plate, and the sensor frame has also been replaced with a wooden version. 
The only eddy currents left in the system axe the currents in the iron cores, and 
these are necessary for the operation of the system. No attempt will be made to 
model the iron core due to the fact that the effect on the system is small, and a 
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limited number of eddy currents can be incorporated in the system using the 
method outlined here. For the purposes of this demonstration the dummy sensor 
ring will be the eddy current that is being modeled. This is because it has the 
most noticeable effect on the system. The eddy current causes the step response 
in pitch to be under damped. 

To show this effect a new controller has been designed, based on a linear 
quadratic regulator controller, that is sensitive to the effects of the eddy current. 
This controller is used for the experimental pitch step responses, and it is also 
used for the simulated pitch step response. By using the same controller an 
amount of uncertainty is eliminated. Figure 4.10 shows a comparison between 
the experimental step in pitch with the dummy sensor ring in place and without 
the ring. By adding the ring the damping of the step response has been greatly 
reduced. 

Figure 4.11 is the same step in pitch, but for this plot the responses were 
simulated using the MATLAB script file located in Appendix C. This script file 
contains a block titled eddy current effects, and that block adds the dynamics of 
the eddy currents to the plant model. The transfer function found to best 
describe the eddy current has been changed to a state-space representation by 
using the MATLAB tf2ss command. The dynamics of this state-space 
representation have been added together, using the MATLAB append command, 
to get a model with five inputs and five output. The five inputs represent the 
five coil currents from the controller, and the five outputs represent the five coil 
currents that have been updated with the eddy current dynamics. The updated 
eddy currents axe fed to the plant model to find the simulated pitch response. 

The change in step response for pitch, between the metal and non metal 
case, shown in Figure 4.11 is not as large a difference as in the case of the 
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experimented step responses, but the trend of making the step response less 
damped with the addition of the dummy sensor ring does appear. The difference 
in the magnitude of the change in damping is not very surprising. The 
theoretical model still has errors, and Figure 4.12 shows the step responses with 
out the dummy sensor ring for the simulation and the experiment. The step 
responses should be exactly the same. Their differences demonstrate the errors 
that remain in the model of the system. One of the known sources of these 
errors are the iron cores which were not added to the dynamics of the plant for 
the purposes of the simulation. 

Gains can be incorporated in the model to make the experiment and 
simulation agree better. The improved simulated step response vs. the 
experimental step response is shown in Figure 4.13, again this plot is created 
without the dummy sensor ring in place. Using this model, which has less initial 
errors than the model that was used originally, the step response was simulated 
for both with and without the dummy sensor ring. This plot is shown in Figure 
4.14, and it agrees better with the experimental measurements shown in Figure 
4.10. 
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Figure 4.5 Dummy Sensor Ring 
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Air Cored Coil and Dummy Sensor Ring 
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Figure 4.7 Response of the Aluminum Plate 
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Figure 4.8 Response of the Iron Core 
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Iron-Cored Coil with Plate and Dummy Sensor Ring 




Figure 4.9 Response of all of the Eddy Currents in the System 


















CHAPTER FIVE 


Magnetic Field Calculation 
5.1 OPERA Calculation Theory 


OPERA is a well established environment for analyzing magnetic 
problems. From this environment a number of different magnetic field solving 
programs, such as GFUN, Elektra, and Soprano, can be run. GFUN was 
developed originally by the Rutherford Appleton Laboratory in England, and its 
current developer and distributor is Vector Fields Inc. it performs static 
magnetic field calculations. It is capable of calculating the effects of a non-linear 
material on the magnetic fields. Most of the theoretical variable calculations did 
not need GFUN for the magnetic field calculation. When the work consisted of 
electromagnets the field calculations could be done in the OPERA post processor 
environment. The only case were GFUN was required is when the iron cores are 
added to the system. 

The fields are calculated in two parts. The first part of the field is from 
the iron cores, and the second part is from the electromagnetic coils as shown 
below. 


H = H- + H 


iron 


1 coils 
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The field from the coils is calculated using Biot Savart law which is integrated 
over the volume of the coil. 


fl = / l# dn ” 

The fields from the iron is not as easy to calculate as the coil’s fields. The iron is 
broken down into small pieces that axe modeled like dipole magnets. The fields 
axe calculated using the equation below: 

H m = V / 

In this equation both H and M axe unknowns, but they are related through the 
equation below. 


M = /i 0 (/*r-l)H 

The previous two equations axe evaluated at each of the control points in the 
finite element grid. This will produce a set of n equations with n unknowns. 

The equations axe solved to find the magnetization of each element. Once the 
strengths of the dipole magnets axe known finding the field at any given point is 
a matter of summing the effects of the dipoles together. 

5.2 Field Gradients 

Calculating the field gradients utilizes the grid command from OPERA. 
The grid command is the only command, in OPERA, that will produce an 
output file that is easily readable by other codes. However, there is no provision 
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in OPERA for finding the field gradients. In order to find the gradients the x, y, 
and z field component must be known along several lines. The configuration of 
these lines is shown in Figure 5.1. Each of these lines has 10 points in it, and 
they are centered around the point of interest. This is the information necessary 
to find all of the field gradients. After the information is retrieved from OPERA 
a mathematical analysis program such as MATLAB is utilized to find the actual 
gradients. For second order gradients, such as B^, a best fit 2 n< ^ order 
polynomial is generated for the plot of B x vs. x. 


B x = A x 2 + B x + C 


The first derivative is 


B xx = 2 A x + B 

evaluated at x=0. For B yx the field values for By along the x axis line as use for 
the best fit polynomial. The second order gradients, such as B xxx , are found 
from the same best fit polynomial. 


Bxxx — 2A 

All of the field gradients can be found from the three on axis lines except for the 
B xyz gradient. The three lines in the x direction are use for this purpose. The 
value of B xy is found on all three lines by the method outlined above. Then a 
2 n< ^ order polynomial is fitted to the three values of B X y giving the B xyz field 
gradient. A listing of the MATLAB script file that is used to find the field 
gradients is listed in appendix B. 
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CHAPTER SIX 


Yaw Angle Interpolation 


Finding the yaw angle of the suspended element with respect to the coils 
is necessary to demonstrate 360 degrees of yaw control. This is due to the fact 
that the computer varies a mixing matrix that determines the amount of current 
each coil should receive in order to maintain suspension. 


I = M * F 


Where I is the coil currents, and F is a vector of the required forces and 
torques. The mixer matrix M changes according to the yaw angle. The problem 
of finding the yaw angle is in the sensors. When the suspended element leaves 
the range of the sensors they do not feedback the correct position, and they have 
a short range. In order to keep the suspended element in sensor range during a 
360 degree rotation the sensors are rotated by hand with respect to the 
electromagnets. This means that the computer needs information as to how far 
the sensors have been rotated. A rotary encoder, or similar piece of equipment, 
could have been used, but a desire not to add any more new hardware to the 
system has led to the development of other methods to detect the yaw angle. 
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Two methods have been developed that use the same amount of information that 
the computer was already receiving. 

In the first method to be implemented, the angle is determined from the 
yaw error [l 7] . As the sensors are rotated, the mixing matrix, which is not being 
updated at this point, develops inaccuracies. These inaccuracies cause the model 
to be suspended with a yaw angle error. When this error is large enough not to 
be noise the controller updates the mixer matrix in such a way as to drive the 
error back to zero. The method has been used to demonstrate 360 degrees of 
yaw, and it is surprisingly reliable. However, the motion is not smooth because 
when the mixer matrix is updated the suspended element “snaps” to the new 
yaw angle. 

The second method reads the time- averaged current going to the coils 
during suspension. For each yaw angle there is a unique current distribution 
required to maintain suspension. The currents in the coils vary like five cosine 
waves shifted in phase as the suspended element is rotated 360 degrees (see 
Figure 6.1). This allows each yaw angle to have a unique set of coil currents. 

By reading the currents an algorithm can determine the angle of yaw and 
evaluate a new mixer matrix. The has the advantages of continuously updating 
the matrix allowing a smooth rotation. 

6.1 The Method of Calculation 

In theory the currents in the coils follow the equation: 

In = Imax * cos(0 n -$) 
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Where I n is the current in the n 1 * 1 coil, I ma x is maximum current 
needed in any coil for the purposes of suspension, 0 n is the angular position of 
the coil itself, and <f> is the position of the suspended element. Both of the angles 
are measured with respect to the first coil. In practice the currents are updated 
at high frequency, to maintain suspension, which causes noise. This noise must 
be eliminated to get an error-free result. The noise is reduced by passing the 
current values through a low pass filter with a break frequency on the order of 
one Hz. 

The yaw function that is generated is as follows: 

5 

Yaw Function = ^ I n * cos(0 n - (tp + e)) 
n=l 

Where ip is the previous yaw angle, and e varies between ± five degrees. 
The yaw function is calculated for each e. When the yaw function is plotted 
versus e, its maximum is at the yaw angle (as shown in Figure 6.2 for a yaw 
angle of 5 degrees and a guessed angles of 3 degrees). 

The distance d is the correction factor to be added to the old value of the 
yaw angle. The distance is calculated be fitting a second order polynomial to the 
yaw function as in the equation below. 

o 

Yaw Function = Ay aw * + By aw * e + Cy aw 

Then the correction factor and yaw angle are: 

■Byaw 

emax - (2 - Ayaw } 


59 



Yaw Angle = 0 + e max 

During initialization of the computer, mixer matrices are stored for six 
degree increments around the first 72 degrees. From these matrices the mixer 
matrix for any angle can be interpolated. Due to the unique symmetry of the 
system the numbers inside the mixer matrix repeat every 72 degrees. The only 
thing that changes is the numbering convention used for the coils. As the 
rotation passes the 72 degree mark coil 2 becomes coil 1 for the purposes of 
calculating the mixer matrix. 

6.2 Experimental Results 


Using the current interpolation program to control the gain scheduling has 
resulted in reliable suspension through a large yaw angle. This program has also 
demonstrated rotations larger than 360 degrees. As the angle passes 360 degrees 
the program resets the angle back to 0 degrees. The rotation using the 
interpolation program for gain scheduling is smoother than when the gain 
scheduling was performed using the yaw angle error filtering routine. Figure 6.3 
shows the suspended element undergoing a large angle rotation. 
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Figure 6.3 Continued 
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CHAPTER SEVEN 


Discussion 

The procedure of talcing experimental frequency response data and using 
it to update parameters in a theoretical formula can be applied to a variety of 
systems such as wind tunnels and magnetic bearings. How accurate the 
mathematical model adjustment is will be determined, in part, by some 
characteristics of the eddy currents involved. When the mathematical model of 
LAMSTF was updated for the aluminum sensor ring the results were very good. 
However, when the model was updated for the iron cores, the results were not as 
accurate, particularly at the higher frequencies. In the case of LAMSTF the 
effect of the cores is not as pronounced as the effects of the sensor ring. 

Therefore, representation of the dominant eddy current leads to significant 
improvements in the model fidelity. 

When applying this procedure to other systems, the eddy current with the 
largest effect must be identified. Once this is done, determination of the useful 
benefit depends on the eddy current itself and the frequencies of interest from 
the controls stand point. Some points to consider about the eddy current itself 
are how well the shape is known or constrained, and if the shape is stable with 
respect to the frequency. 

In the case of the iron core the shape of the eddy current changes as the 
frequency of the driving magnetic field is increased. This is due to the change of 
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the skin depth, which has otherwise been assumed to be so large that the eddy 
currents are distributed evenly throughout the piece of metal. As the shape 
changes, the resistance, which was assumed to be constant, is drastically altered. 
Checking the skin depth requirement at high frequency will give an indication of 
the applicability of the technique to other eddy currents. 

Finding the shape of the eddy current is important to establish variables 
that can be used as a starting point for estimating the parameters that will then 
be used experimentally to update the system’s mathematical model. The more 
that is known about the shape of the eddy current the more accurate the starting 
point will be. If the eddy current is highly constrained, as was the case of the 
sensor ring, the initial guess will be fairly close to the actual answer. In the case 
of the aluminum plate, the eddy current’s diameter was unknown, and this made 
the first approximation worse as compared to the first approximation of the the 
sensor ring. The same thought holds true for multiple eddy currents. If the 
eddy currents can be isolated, and then added together the theory will work. If, 
however, the eddy current can not be isolated the effects of a single eddy current 
will be lost when combined with all of the other eddy currents. This was the 
case with the original sensor support frame where one piece of metal could be 
part of 3 or 4 different eddy current loops. 

The last point to consider is the frequency range of interest. The 
mathematical model does not have to be accurate out to a high frequency if it is 
not used. The range for frequencies can be determined by the highest open loop 
unstable pole of the system. This pole will determine how fast a controller must 
be run to achieve suspension. If the analysis corresponds to the experimentally 
determined transfer function then the updates to the model will be helpful. 
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CHAPTER EIGHT 

Conclusions 


When the system was under construction, metal was used extensively for 
support structures, etc.. It was known that eddy currents would be present in 
the system, and the purpose of using so much metal was to study the effects of 
these eddy currents. The eddy currents did not effect the system to the point 
where the system could not be made operational. They did, however, have a 
pronounced effect on the performance of the system. This effect showed up in 
the step responses of the different degrees of freedom most notably in the pitch 
response. After the general effect was understood, the system was refined to 
study the individual eddy currents. This was done by eliminating the eddy 
currents using traditional techniques such as replacing the metal with a 
nonconductive material. This step was necessary to isolate the different eddy 
currents. Before the system was altered it was not clear which eddy currents 
were altering the system dynamics. After the different eddy currents were sorted 
out the sensor frame was shown to have the most prominent effect on the system 
followed by the iron cores, and the aluminum baseplate effect was shown to be 
very small. A means to add the dynamics of the sensor ring to the system was 
developed, and it is outlined in this thesis. The theory best applies to a dummy 
sensor ring. That case was used to compare the updated model to the 
experimental system. By adding the ring’s dynamics to the system the step 
response in pitch became less damped, and this was the effect observed in the 


experiments on the system. 
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The method outlined for analyzing the eddy currents to the system has 
several drawbacks. One eddy current, the dummy sensor ring, has been added to 
the system. The results agreed better with the experiment, but errors still exist 
in part due to the other eddy currents. These other eddy currents were not 
added to the system for two reasons. The theory did not cover the iron cores as 
well as the sensor ring, and the method by which the eddy currents are added is 
questionable when a significant amount of eddy current need to be added to the 
system. Adding together several different eddy currents is necessary to improve 
a system with a lot of metal and multiple eddy currents in the structure. A 
better method is needed for adding multiple eddy currents to the system 
dynamics. 

Another draw back to the theory is the problem of describing the eddy 
currents present in the iron cores. A theory is under development that can 
explain the dynamics of the iron cores. One of the unique problems of describing 
the eddy currents present in the core is that the gain roll off is 10 db/dec, or half 
of a single pole’s rolloff. This is explained in the proposed theory by using an 
infinite pole zero arrangement. This will be difficult to add to the system due to 
the fact that even when truncating to just 5 or 6 poles the method currently 
being used to add the dynamics to the system becomes questionable. 
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APPENDIX A 


Code Development 

The code was developed based on a code that was developed for the 
original system. The original code is discussed in more detail in reference 17. 
Severed modifications were performed on the original code due to updates in the 
system. These updates include changing the boards used, deeming up the code 
and installing the yaw angle interpolation subroutines. 

Updating the program is necessary when a different board is used because 
every board heis different commands that serve the different input-output ports. 
To help facilitate the updating and correction of the program it has been 
moduralized. This procedure starts with the original program that was 26 pages 
of code and breaks it into several different files each of which control a single 
aspect of the program’s performance. This makes the program easier to update 
leaving the bulk of the program intact. A flow chart for when the different files 
are used is presented followed by a listing of the program. 
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Header file to declare vaxiables and subroutines. 


^define TIMER0x210 /* Timer board Address */ 

^define DD Abase 0x320 /* D/A base address */ 

#define UP 72 

#define DOWN 80 

#define LEFT 75 

# define RIGHT 77 

#define ENTER 28 

#define ESC 27 

^define FI 59 

#define F2 60 

#define F4 62 

^define F5 63 

^define F6 64 

^define F9 67 

#define F10 68 

^define pi 3.14159 

^define On 1 

^define Off 0 

#define MAXCUR 8.36 

^pragma intrinsic(outp,inp) 

#i£ndef WHERE 

^define WHERE extern 

WHERE double theta[5]; 

WHERE unsigned int chan; 

WHERE int rota; 

WHERE int chans[5]; 

WHERE int intnum; 

WHERE int gains [5]; 

WHERE int plotcolor,datcolor; 

WHERE int zoom,ang; 

WHERE int tmp[5],ix; 

WHERE double offset [5]; 

WHERE double gain[5],gainstep,positionstep; 

WHERE double mix[5][5]; 



WHERE double yawangle,yawanglestep, autoyaw; 

WHERE int chd; 

WHERE char plotmode, ch2; 

#else 

double theta[5]={0.0, 72.0, 144.0, 216.0, 288.0}; 
int chd=l; 

char plotmode= , a’,ch2=’z’; 
unsigned int chan=10; 
int rota=0; 

int chans[5]={0,l,2,3,4}; 

int intnum=0x0d; /* interrupt vector number */ 
int gains[5]={0,0,0,0,0}; 
int plotcolor=0,datcolor=15; 
int zoom=l,ang=0; 

int tmp [5]={53,101,149,197,245},ix=25; 
double o£fset[5]={8.36, 2.713,-6.965, -6.965, 2.713}; 
double gain [5] ={1,1,1, .5, .5} ,gainstep =.01 ,positionstep= . 
double mix[5][5]={{ 1.0, 0.0, -1.0, 0.625, 0.0}, 
{-0.809,0.618,-0.309,1.0,1.0}, 

{ 0.309,-1.0,0.809,0.768,0.618}, 

{ 0.309,1.0,0.809,0.768,-0.618}, 
{-0.809,-0.618,-0.309,1.0,-1.0}}; 
double yawangle= 0 ,y awanglestep= 0 . 5 , autoyaw=0; 
#endif 

WHERE int mode, flag, params[10]; 

WHERE unsigned int *DMA_buf,*bufoffset; 

WHERE int far *chans_ptr; 

WHERE int far *gains_ptr; 

WHERE double mixail [5] [5] [6 1 ] ,filt [5] , positional; 
WHERE unsigned int adcsr; 

WHERE int kv, dec, sampler; 

WHERE int yll,xl,samplrate,diviser; 

WHERE int input, k,j,num; 

WHERE int ii,i,l,x,olddat [500] [5]; 

WHERE int positionsign,gainupdown,trig; 

WHERE double all,a,b,c,d,e,f,current[5],savdat[5][500); 



WHERE double v,vl[25],v2[25],vsensor[6]; 
WHERE double cosdat[800]; 

WHERE int maxnum,yawerror limit; 

WHERE chax ch, chi, loop, chin; 

WHERE long backc; 

WHERE double y,xy,x2y; 

WHERE double correction, testangle, function, ag; 
void (_interrupt _far *oldnum)(void); 
void _interrupt _far isr(void); 
extern MSCL_DAS40(int*,int*,int*); 

void initialize(void); 
void oscilloscope(void); 
void timerset(void); 
void adcinitialize(void); 
void plotdat(void); 
void keyboard(void); 
void outdat(void); 
void gainplot(void); 
void positionplot(void); 
void compensator(void); 
void timersetup(void); 
void acknowledge (void); 
void savedat(void); 
void getdat(void); 
void newmixer(void); 
void cosmak(void); 
void sbutdown(void); 
void isrsetup(void); 
void yawer(void); 



Interrupt Service routine (controller) 


^include <graph.h> 

^include <stdio.h> 

^include <conio.li> 

^include <math.h> 

^include <float.h> 

^include <dos.h> 

^include <bios.h> 

^include <signal.h> 

^include <stdlib.h> 

^include <malloc.h> 

^include ” suspend. h” 

/* Intrrupt service routine */ 

/* datin() obtains data from ADC in the array */ 
void ..interrupt _far isr(void) 

{ 

unsigned int yll,Outport; 
int k,df,tmp; 
for(k=0;k<=5;k+-f) 

{ 

mode = 7; 
flag = 0; 
params[0] = k; 
params[l] = 0; 

MSCL_DAS40 (&mode, params, &flag); 
if (flag != 0) ch=ESC; 
tmp=params[0]; 
vsensor[k]=(float)(tmp/409.5); 

} 

_enable(); 

I ********** ************************************************************ 

*/ 

/* decoupler */ 



vl[2]=((vsensor[2]-vsensor[3])+(vsensor[5]-vsensor[4]))/(4*.707); 
vl [5] = ((-vsensor[2]+ vsensor [3] )+ (vsensor [5]- vsensor [4]))/(4*. 707); 
v 1 [8] = ( vsensor [0] + vsensor [ 1 ] ) / 2 ; 
vl[H]=( vsensor [ 1 ] - vsensor [0] ) ; 

vl[14]=((-vsensor[2]-vsensor[3])+(vsensor[4]+vsensor[5]))/2; 

if(sampler<=200) 

{ 

for(df=0;df<=4;df++) 

savdat [df] [sampler] = vl [2+df*3] ; 

sampler++; 

} 

/* compensator: dual phase advance using backward difference method */ 

/* vl[2+3*l]=vl[2+3*l]+position[l]*4; */ 

for(l=0;l<=4;l++) 

{ 

v2[2+3*l]=all*(gain[l]*(b*(vl[2+3*l]+position[l]) 

+c*vl[l+3*l]-fd*vl[3*l])- 

e*v2[l+3*l]-f*v2[3*l]); 

} 

I*********************************************************************** 

*1 

/* current formation*/ 
for(l=0;l<=4;l++) 

{ 

current [I]=v2[2]*mix[l][0]+v2[5]*mix[l][l]+v2[8]*mix[l][2]+ 
v2 [ 1 1 ] * mix [1] [3] + v2 [ 1 4] *mix [1] [4] +off set [1] ; 

/******************************************************************/ 

/* Data outputs to the five D/As */ 

if (f abs ( current [1] ) < 1 0 . 0 ) 
yll = (unsigned int)(204.75*(current[l]+10.0)); 

Outport=DDAbase-(-l*2; 

_disable(); 

_asm 

{ 

push ax; 
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push dx; 
mov dx,Outport; 
mov ax,yll; 
out dx,ax; 
pop dx; 
pop ax; 

} 

_enable(); 

vl[3*l]=vl[l+3*l]; 
vl[l+3*l]=vl [2+3*1]; 
v2[3*l]=v2[l+3*l]; 
v2[l+3*l]=v2[2+3*l]; 

} 

_disable(); 

_asm 


{ 

push ax 
push dx 

mov dx, TIMER ;Load the ADC board address 

add dx,6 ;point to adc high byte register 
or al,al ; clear al 
in al,dx ;output command 
mov dx,20h ; load 8259 address 
mov al,20h ; set normal priority, EOI=l, (65h) 
out dx,al ; output command 
pop dx 
pop ax 


_enable(); 
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Main Program 


^define WHERE 
^include <graph.h> 

^include <stdio.h> 

^include <conio.h> 

^include <math.h> 

^include <float.h> 

^include <dos.h> 

^include <bios.h> 

^include <signal.h> 

^include <stdlib.h> 

^include <malloc.h> 

^include ” suspend. h” 

void main(void) 

{ 

/* Start up and initilization routines */ 

isrsetup();/* sets up the interrupt vector */ 
cosmak();/* calculate cosine values for every .5 deg */ 
getdat();/* Get mixer matrix data */ 
adcinitialize( ) ; /* Initialize the A/D converter board */ 
initialize();/* Compute The compensators parameters */ 
oscilloscopeQ;/* Display the Controller Screen */ 
timersetup();/* Set timer to interrupt at sample- rate */ 

/* This is the controller loop */ 

ch=ESC; 

while(chl!=ESC) 

{ 

_setcolor(4); 

_rectangle(_GFILLINTERIOR,550,420,637,477); 

_settextposition(28,71); 

P rintf(” RUNNING”); 

_settextposition(29,71); 
printf(”ESC- quit”); 



outp(TIMER+l,0x61);/*Starts the interrupt by loading and arming timers*/ 
while(ch!=ESC) 

{ 

keyboard(); /* Check for user input */ 
yawer(); /* This does the gain scheduling */ 
for(j=0;j<=4;++j) 

{ 

v=vl[2+3*j]; 

if(sampler<=2) 

{ 

_setcolor(14); 

_rectangle(_GFILLINTERIOR,600,330,620,325); 

} 

else 

if(sampler>200) 

{ 

_setcolor(l); 

_rectangle(_GFILLINTERIOR,600,330,620,325); 

} 

if(chd!=0) 

{ 

if(!trig || ix<474) 

{ 

if(plotmode==’a’) plotdat(); 
if(plotmode==V) 

{ 

switch(ch2) 

{ 

case ’x’: if(j==0) plotdat(); 
break; 

case ’y’: if(j==l) plotdat(); 
break; 

case ’z’: if(j==2) plotdat(); 
break; 

case ’m’: if(j==3) plotdat(); 
break; 



case ’n’: if(j==4) plotdat(); 
break; 

} 

} 

} 

} 

} 

if((trig==0)&&(ix>=475)) ix=25;/* puts oscope counter back at start */ 
else ix-f-f; /* or advances it */ 

_setcolor(2); 

_rectangle(_GFILLINTERIOR,550,420,637,477); 

_settextposition(28,71); 

_outtext(”R-Run”); 

_settextposition(29,71); 

_outtext(”ESC-exit”); 

outp(TIMER+l,Oxdf); /* shut down interrupt counter */ 
for(ii=0;ii<=4;++ii) 

{ 

outp(DDAbase+ii*2,0); /* 0000000 for low byte at zero */ 
outp(DDAbase+ii*2+l,8); /* 1000 for high byte at zero */ 

} 

ch=V; 

chl=(char)getch(); 

} 

_setvideomode(_DEFAULTMODE); 

_dos_setvect(intnum,oldnum); / * reset interrupt vector table */ 
system^ d40 u”); 



Interrupt Setup 


^include 
^include 
# include 
^include 
^include 
^include 
^include 
^include 
^include 
^include 
^include 


<graph.h> 

<stdio.h> 

<conio.h> 

<math.h> 

<float.h> 

<dos.h> 

<bios.h> 

<signal.h> 

<stdlib.h> 

<maIloc.h> 

”suspend.h” 


void isrsetup(void) 

{ 

oldnum=_dos_getvect(intnum); /* OxOd IRQ-5 */ 

_disable(); /* Disable all interrupts */ 

_dos_setvect(intnum,isr); / * Set new interrupt vector table */ 

_enable(); /* Enable all interrupts */ 
outp(0x21,0x00); /* Send acknowledgement signal */ 

outp(0x20,0x20); /* To the interrrupt handler */ 


} 



Cos Calculator 


^include <graph.h> 
^include <stdio.h> 

# include <conio.h> 

# include <math.h> 
^include <float.h> 
^include <dos.h> 
^include <bios.h> 
^include <signal.h> 
^include <stdlib.h> 
^include <malloc.h> 
^include ” suspend. h” 


/* cosmak : calculates cos values for angles 0..360 */ 

void cosmak(void) 

{ int i; 

double radang; 

maxnum= (int ) (360 . 0/ yawanglestep) ; 
yawangle=0; 

for(i=0;i<=maxnum;i++) 

{ 

radang= (yawangle / 1 80 . 0) *pi ; 
cosdat [i] =cos(radang) ; 
yawangle=yawangle+yawanglestep; 

} 

yawangle=0; 



Retrives the Mixer Matrix Data 


^include <graph.h> 

^include <stdio.h> 

^include <conio.h> 

^include <math.h> 

^include <float.h> 

^include <dos.h> 

^include <bios.h> 

^include <signal.h> 

^include <stdlib.h> 

^include <ma!loc.h> 

^include ” suspend. h” 

/* getdat() loads mixing matrix data*/ 
void getdat(void) 

{ 

int ij,k,dim; 

FILE *fp; 

printf(” Getting Data From MIX.DAT\n”); 
if((fp=fopen(”mix.dat”,”r+”))!=NULL) 

{ fscanf (fp,” %d” ,&dim); 
printf(” Number of matrices axe: %d \n”,dim); 
for(k=0;k<=dim-l;k++) 

{ 

for(i=0;i<=4;i++) for(j=0;j<=4;j++) fscanf(fp,”%lT,&mixall[i][j][k]); 

} 

printf(”Data was successfully loaded\n”); 
fclose(fp); 

} 

else 

fprintf(stderr,” Oops file error”); 

} 
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ADC board initilizer 


^include 
^include 
^include 
^include 
# include 
^include 
^include 
^include 
^include 
^include 
^include 


<graph.h> 

<stdio.h> 

<conio.h> 

<math.h> 

<float.h> 

<dos.h> 

<bios.h> 

<signal.h> 

<stdlib.h> 

<malloc.h> 

"suspend.h” 



void adcinitialize(void) 

{ 

/* DAS-40 Initialization */ 

system(”d40”); 

mode=0; 


flag=0; 

params[0]=0; 

MSCL_DAS40(&mode,params,&flag); 
if(flag !=0) 

{ 

printf (”**** Error %u detected in mode %u ”, flag & Oxff, ((flag & OxffOO) 
» 8 )); 
exit(l); 

} 


} 
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Initilizes the controller parameters 


^include <graph.h> 
^include <stdio.h> 
^include <conio.h> 
^include <math.h> 
^include <float.h> 
^include <dos.h> 
^include <bios.h> 
^include <signal.h> 
^include <stdlib.h> 
^include <malloc.h> 
^include ” suspend. h” 


/* initialize: setting the parameters of the dual phase advance compensator */ 
void initialize(void) 

{ 

double dt,T,n; 
int breakfreq; 

printf(”\nSample rate : 500\n”); 

samplrate=500; /* Get the sampling rate */ 

printf(” Compensator’s Break frequency is 170hz\n”); 

breakfreq=170; /* Get the break frequency */ 

dt=l /((double) samplrate); /* period of the sampling */ 

T=l/(2.0*pi* (double) (breakfreq)); /* 1/ break frequency*/ 
printf(”\nT= %9.6f\n\n”,T); 

n=6.0; /* ratio of lead to lag break frequency */ 

all=l/((dt*dt)+4*T*dt+4*(T*T)); 

b=(dt*dt)+4*n*dt*T+4*n*n*T*T; 

c=2*dt*dt-8*n*n*T*T; 

d=(dt*dt)-4*n*dt*T+4*n*n*T*T; 

e=2*dt*dt-8*T*T; 

f=(dt*dt)-4*dt*T+4*T*T; 

printf(”The compensator parameters are:\n”); 



printf(”a=%9.6f\n n ,all); 
printf(”b=%9.6f\n” ,b); 
printf(”c=%9.6f\n ,, ,c); 
printf(” d=%9.6f\n” ,d) ; 
printf(”e=%9.6f\n”,e); 
printf ( ” f = %9 . 6f\n” ,f ) ; 

getch(); 



Draw the Screen Display 


^include <graph.h> 

^include <stdio.h> 

^include <conio.h> 

^include <math.h> 

^include <float.h> 

^include <dos.h> 

^include <bios.h> 

^include <signal.h> 

^include <stdlib.h> 

^include <malloc.h> 

^include ” suspend. h” 

/ * oscilloscope: signal display screen setup */ 
void oscilloscope(void) 

{ 

int lljj,gainstart,posibar; 

_ set videomode( _ VRES 16C0L0R); 

/* initialization of the screen clear array */ 
for(ll=0;ll<=5;++ll) 

{ 

for(jj=0;jj<=500;++jj) 

olddat[jj][ll]=20; 

} 

I* display of the oscilloscopescreen */ 
_clearscreen( _ GCLEARSCREEN ) ; 
_setcolor(13); 

_rectangle(_GBORDER,0,0,639,479); 

_setcolor(plotcolor); 

_rectangle(_GFILLINTERIOR,20,5,480,290); 

_settextposition(4,2); 

printf(”X”); 

_settextposition(7,2); 

printf(”Y”); 

_settextposition(10,2); 



printf(”Z”); 

_settextposition(13,2); 

printf(”M”); 

_ set textposi tion ( 1 6 ,2 ) ; 
printf(”N”); 

_setcolor(7); 

for(U=52;ll<=244;ll=ll+48) 

{ 

_moveto(19,ll); 

_lineto(24,ll); 

_moveto(476,U); 

_lineto(484,U); 

} 

/*^**** ********************************♦**************/ 
/* Gain control knobs Display */ 

_setcolor(12); 

_rectangle(_GBORDER,2,300,480,478); 

_settextposition(20,2); 

_settextcolor(14); 

_outtext(”Feedback gain control”); 

_settextposition(20,40); 

_outtext(” Position Control”); 

_setcolor(5); 

jj=0; 

for(ll=21;ll<=29;ll=ll+2) /* lableing the gain knobs */ 

{ 

_settextposition(ll,27); 
printf (” %4. 3F ,gain [j j ] ) ; 

jj++; 

_settextposition(ll,4); 

printf(”0”); 

} 

for(ll=319;ll<=447;ll=ll+32) j* gain control bar display */ 

_ rect angle ( _ G B 0 RD ER,42 ,11 , 1 96 ,11 +22 ) ; 



jj=0; 

for(U=320;ll<=448;ll=ll+32) /* Presetting gain display */ 

{ /* on the screen */ 

gainstart=(int)(gain[jj]*80)+40; 

_setcolor(6); 

_rectangle(_GFILLINTERIOR,40,ll,gainstart+3, 11+20); 

++jj; 

} 

/* display of position control knobs */ 


_setcolor(4); 

for(U=320;ll<=448;ll=ll+32) /* position control bar display */ 

_rectangle(_GBORDER, 310, 11, 460, 11+20); 

_setcolor(10); 

jj=0; 

for(U=320;ll<=448;ll=ll+32) 

{posibax=385+(int)(position[jj]*75); 

_rectangle(_GFILLINTERIOR,posibar,ll,posibar+l,ll+20); 


++jj; 

} 


/***************************************/ 


/* labling of gain knob display */ 


_settextposition(21,2); 

_outtext(”X”); /*x*r axial*/ 
_settextposition(23,2); 

_outtext(”Y”); /* y or side */ 
_settextposition(25,2); 

_ out text (”Z”); /* z or vertical */ 

_ settextposi tion ( 2 7 ,2 ) ; 

_outtext(”M”); /* M - pitch movement */ 
_settextposition(29,2); 

_outtext(”N”); /*N - yaw movement */ 

y ************************************/ 


/* Display of options */ 
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_settextposition(2,62); 
printf(”E- enable plot”); 
_settextposition(3,62); 
printf(”D- disable plot”); 
_settextposition(5,62); 
printf(” Select Channel”); 
_settextposition(6,62); 
printf(” Press ”); 
_settextposition(7,62); 
printf(” X Y Z M N”); 
_settextposition(9,62); 
printf(”S- single plot”) ; 
_settextposition( 10,62); 
printf(”A - all plots”); 
_settextposition( 1 1 ,62) ; 
printf(”T - triggered ”); 
_settextposition(12,62); 
printf(”C - continuous”); 
_settextposition(13,62); 
printf(”Fl- unzoom plot”); 
_settextposition(14,62); 
printf(”F2- zoom plot”); 
_settextposition( 15,62); 
printf(”F4- Save data”); 
_settextposition( 16,62); 
printf(”F5- - Yaw ”); 
_settextposition( 17,62); 
printf(”F6- + Yaw ”); 

_ settextposi tion ( 1 8 ,62 ) ; 
printf(”F9- AutoYaw off”); 
_settextcolor(4); 
_settextposition(18,75); 
_outtext ( ” off” ) ; 

_ settextposi tion (19,62); 
printf(”F10- AutoYaw on”); 
_settextposition(24,62); 



printf(”Yaw Error:”); 
_settextposition( 22 , 62 ); 
printf(”Yaw Angle:”); 

} 



Timer board setup 


# include 
^include 
^include 
^include 
^include 
^include 
^include 
^include 
^include 
^include 
^include 


<graph.h> 

<stdio.h> 

<conio.h> 

<math.h> 

<float.h> 

<dos.h> 

<bios.h> 

< signal. h> 
<stdlib.h> 
<malloc.h> 
” suspend. h” 


/* timersetup: sets up the sample rate clock on the second ADC for interrupt */ 
void timersetup(void) 

{ 

int rate; 

/* Timer ic AM9513 setting */ 

outp(TIMER+l,Oxdf); /* disarms counters */ 
outp(TIMER+l,23); /* sets next outp to the master register*/ 
outp(TIMER,OxcO); /* master rgister low byte */ 
outp(TIMER,0x41); /* master register high byte */ 
outp(TIMER+l,0x01); /* sets next outp to the counter 1 register */ 
outp(TIMER,0x22); /* sets wave form and base frequency */ 
outp(TIMER,0x0c) ; 
rate=(int)(31250/samplrate); 

outp(TIMER+l,0x09); /* sets next outp to the load register */ 
outp(TIMER,rate); /* it will count down from rate and set off an interrupt 

*/ 

outp(TIMER+l,0x61); /* load and arm */ 


97 



Tests for Keyboard Interaction 


^include <graph.h> 

^include <stdio.h> 

^include <conio.h> 

^include <math.h> 

^include <float.h> 

^include <dos.h> 

^include <bios.h> 

^include <signal.h> 

^include <stdlib.h> 

^include <maIloc.h> 

^include ”suspend.h” 

/* keyboard: User interaction through keyboard*/ 
void keyboard(void) 

{ 

if(kbhit()!=0) 

{ 

chin=(char)getch(); 

switch(chin) 

{ 

case ’X’rcase ’x’: ch2=’x’; 

_setcolor(8); 

_rectangle(_GFILLINTERIOR,35,320,39,478); 

_setcolor(12); 

_rectangle(_GFILLINTERIOR,35,320,39,340); 
if(plotmode= = ’s’ ) 

{ 

_ setcolor (plotcolor) ; 

_rectangle(_GFILLINTERIOR,20,5,480,290); 

} 

break; 

case ’Y’rcase ’y’: ch2=’y’; 

_ set color (8); 

_rectangle(_GFILLINTERIOR,35,320,39,478); 



_setcolor(12); 

_rectangle(_GFILLINTERIOR,35,352,39,372); 

if(plotmode==’s’) 

{ 

_setcolor(plotcolor); 

_rectaiigle(_GFILLINTERIOR,20,5,480,290); 

} 

break; 

case ’Z’rcase V: ch2=’z’; 

_setcolor(8); 

_rectangle(_GFILLINTERIOR,35,320,39,478); 

_setcolor(12); 

_rectangle(_GFILLINTERIOR,35,384,39,404); 

if(plotmode==’s’) 

{ 

_ setcolor (plot color ) ; 

_rectangle( _ GFILLINTERIOR,20,5,480,290) ; 

} 

break; 

case ’M’:case ’m’: ch2=’m’; 

_setcolor(8); 

_rectangle(_GFILLINTERIOR,35,320,39,478); 

_setcolor(12); 

_rectangle(_GFILLINTERIOR,35,416,39,436); 

if(plotmode==V) 

{ 

_setcolor(plotcolor); 

_rectangle(_GFILLINTERIOR,20,5,480,290); 

} 

break; 

case ’N’xase ’n’: ch2= , n , ; 

_setcolor(8); 

_rectangle(_GFILLINTERIOR,35,320,39,478); 

_setcolor(12); 

_rectangle(_GFILLINTERIOR,35,448,39,468); 
if ( plotmode= = ’s ’ ) 



{ 

_ set color ( plotcolor ) ; 

_rectangle(_GFILLINTERIOR,20,5,480,290); 

} 

break; 

case ’S’: case ’s’: plotmode=V; 
_setcolor(plotcolor); 

_rectangle(_GFILLINTERIOR,20,5,480,290); 

break; 

case ’A’: case ’a’: plotmode=’a’; 
break; 

case ’T’: case ’t’:trig=l; 
break; 

case ’C’: case ’c’: trig=0; 
break; 
case 0: 

chin=(char)getch(); 

switch(chin) 

{ 

case RIGHT: gainupdown=l; 
gainplot(); 
break; 

case LEFT: gainupdown=-l; 
gainplot(); 
break; 

case UP: positionsign=l; 
positionplot(); 
break; 

case DOWN: positionsign=-l; 
positionplot(); 
break; 

case FI: if(zoom>l) zoom = zoom- 1; 
_settextposition( 1,60); 
printf(”%i ”,zoom); 
break; 

case F2: if(zoom<10) zoom=zoom+l; 



_settextposition(l,60); 
printf(”%i ”,zoom); 
break; 

case F4:savedat(); /* Save step responses if needed */ 
break; 

case F5:yawangle=yawangle-yawanglestep; 
newmixer(); 
break; 

case F6:yawangle=yawangle+yawanglestep; 
newmixerQ; 
break; 

case F9: rota=0; 

_ settextcolor ( 8) ; 

_settextposition(19,75); 

_outtext(”on”); 

_settextcolor(4); 

_settextposition( 1 8 , 75) ; 

_outtext(”off’ ); 
break; 

case F10: rota=l; 

_ settextcolor ( 8 ) ; 

_settextposition(18,75); 

_ out text (” off 1 ); 

_ set text color ( 4 ) ; 

_settextposition(19,75); 

_outtext(”on”); 

break; 

} 

break; 

case ’E’:case ’e’: chd=l; 
break; 

case ’D’xase ’d’: chd=0; 
break; 

case ESC: ch=ESC; break; 

} } 

chin=’.’;} 
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Yaw Interpolation Program 


^include <graph.h> 
^include <stdio.h> 
^include <conio.h> 
^include <math.h> 
^include <float.h> 
^include <dos.h> 
^include <bios.h> 
^include <signal.h> 
^include <stdlib.h> 
^include <malloc.h> 
^include ” suspend. h” 


void yawer(void) 

{ 

int loopl,loop2; 


y=0.0; 

xy=0.0; 

x2y=0.0; 


for(loopl=0;loopl<=10;loopl++) 

{ 

function=0.0; 

for(loop2=0;loop2<5;loop2+-f-) 

{ 

_ set textposition ( 25,0 ) ; 

ag=(double)((theta[loop2]-testangle+(double)(loopl-5)) 
*(double)3.14/(double) 180.0); 

function= ( current [loop2] *cos (ag) ) + function; 

} 



y=y+function; 

xy=xy+function*(loopl-5); 

x2y=x2y+(double)((loopl-5)*(loopl-5))*function; 

} 


correction=((42.9*xy)/(110.0*y-11.0*x2y)); 
if(fabs(correction)>5) correction=0.0; 

ag=testangle - correction; 

testangle=.97*testangle + .03*ag; 

if((testangle >=360)||(testangle<=-360)) testangle=0; 

_settextposition(24,72); 
printf(”%4.11f ”,testangle); 
autoyaw=testangle-yawangle; 

if(rota==l) /* Check for automatic yaw sensing */ 

{ /* option */ 
if(autoyaw>=.5) 

{ 

yawangle=yawangle+yawanglestep; 

newmixer(); 

} 

if ( autoyaw <=-. 5 ) 

{ 

yawangle=yawangle-yawanglestep; 

newmixer(); 

} 

} 
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Calculates the New Mixer Matrix 


^include 

^include 

^include 

^include 

^include 

^include 

^include 

^include 

^include 

^include 

^include 


<graph.h> 

<stdio.h> 

<conio.h> 

<math.h> 

<float.h> 

<dos.h> 

<bios.h> 

<signal.h> 

<stdlib.h> 

<malloc.h> 

’ , suspend.h” 


/* newmixer : calculates new mixing matrix and current distribution */ 
void newmixer(void) 

{ int i,j,nextone,indexl,index2,index3,index4,index5,seventy2; 
int angindex; 

double angstep=6, yawangletemp; 
double interpol; 

angindex=(int)(yawangle/yawanglestep);/* Index of the yaw angle */ 
seventy2=(int)(72.0/yawanglestep);/* Index of 72 deg based on*/ 
nextone=l;/* the yaw step chosen */ 
if (angindex<0) 

{ 

angindex=maxnum+angindex; 

nextone=-l; 

} 

indexl=angindex; /* determination of cosine matrix */ 
index2=abs(seventy2-angindex); /*index for calculating the zero */ 



index3=abs(2*seventy2-angindex); /*position current. Every index* */ 
index4=abs(3*seventy2-angindex); /*is 72 degrees apart from its */ 
index5=abs(4*seventy2-angindex); /*adjasent neighbor. */ 

if(indexl>maxnum) /* The five if statements rap around */ 
indexl=indexl-maxnum; 

if(index2>=maxnum) /* the indceis to the begining after */ 
index2=index2-maxnum; 

if(index3>=maxnum) /* 360 degrees rotation */ 
index3=index3-maxnum; 
if(index4>=maxnum) 
index4=index4-maxnum; 
if(index5>=maxnum) 
index5=index5-maxnum; 


offset[0]=MAXCUR*cosdat[indexl]; 
offset [1 ] =MAXCUR*cosdat [index2] ; 
offset [2] =MAXCUR*cosdat[index3]; 
offset [3] =MAXCUR*cosdat[index4]; 
offset [4] =M AXCUR* cosdat [index5] ; 


/* Levitation current calculation */ 
/* using the cosine array, which */ 

/* selected depending on the angle */ 
/* of the model.lnitialoffset[0] */ 

/* is the highest steady current. */ 


yawangletemp=fabs(yawangle); 


if(yaw«mgle>=360) /* Adjustment for angles>360*/ 
yawangletemp=yawangle-360; 

ang=(int)(yawangletemp/angstep);/* Determination of number of steps */ 
nextone=l; 

/* Determination of the interpolating coefficent */ 

interpol= ( (y awangletemp )- (double) ( ( ang) *angstep) ) / angstep ; 


if(yawangle<0) 

{ 

nextone=-l; 

ang=60-ang; 
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} 


/* Interpolation between the mixing matrices */ 
for(i=0;i<=4;i++) 

{ 

for(j=0;j<=4;j++) 

{ 

mix[i] [j] =interpol* (mixall[i] [j] [ang+nextone] 
-mixall[i][j][ang]) 

+mixall[i][j][ang]; 

} 

} 

_settextposition(22,72); 

printf(” %5.2P ,yawangle) ; 


} 



Graphical presentation of the degrees of freedom 


^include < graph. h> 

^include <stdio.h> 

^include <conio.h> 

^include <math.h> 

^include <float.h> 

^include <dos.h> 

^include Chios. h> 

^include <signal.h> 

^include <stdlib.h> 

^include <malloc.h> 

^include ” suspend. h” 

/* gainplot: to display gain on each degree of freedom graphically */ 
void gainplot(void) 

{ 

int gainknob,knobcolor; 
switch(ch2) 

{ 

case ’x’: 

if((gain[0]>=0) kk (gain[0]<=2)) 

{ 

knobcolor=3* (gainupdown-t- 1 ); 

_ set color ( knobcolor ) ; 
gainknob=40+ (int ) ( gain[0] *80 . 0) ; 

_rectangle(_GFILLINTERIOR,gainknob,320,gainknob+3,340); 

gain[0]=gain[0]+gainupdown*gainstep; 

_settextposition(21 ,27 ) ; 
printf(” %4.3r ,gain[0] ); 

} 

else 

if(gain[0]<-5) gain[0]=-5; 
if(gain[0]>5) gain[0]=5; 
break; 



case ’y’ : 

if((gain[l]>=0) && (gain[l]<=2)) 

{ 

_setcolor(knobcolor); 

gainknob=40+(int)(gain[l]*80.0); 

_rectangle(_GFILLINTERIOR,gainknob,352,gainknob+3,372); 
gain [ 1 ] =gain [ 1 ] +gainupdown *gainstep ; 

_ settext position (23,27); 
printf(” %4.3F ,gain[l] ) ; 

} 

else 

if(gain[l]<0) gain[l]=0; 
if(gain[l]>2) gain[l]=2; 
break; 
case V: 

if((gain[2]>=0) && (gain[2]<=2)) 

{ 

knobcolor=3*(gainupdown+l); 

_setcolor(knobcolor); 

gainknob=40+(int)(gain[2]*80.0); 

_rectangle(_GFILLINTERIOR,gainknob,384,gainknob+3,404); 
gain [2] =gain [2] +gainupdown*gainstep; 

_ settext position (25 ,27) ; 
printf ( ” %4 . 3f” ,gain [2] ) ; 

} 

else 

if(gain[2]<0) gain[2]=0; 
if(gain[2]>2) gain[2]=2; 
break; 
case ’m’: 

if((gain[3]>=0) && (gain[3]<=2)) 

{ 

knobcolor=3 * (gainupdown+ 1 ) ; 

_setcolor (knobcolor) ; 
gainknob=40+(int)(gain[3]*80.0); 

_rectangle(_GFILLINTERIOR,gainknob,416,gainknob+3,436); 
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gain [3] =gain [3] +gainupdown*gainstep; 

_settextposition(27,27); 
printf(” %4.3f” ,gain[3] ) ; 

} 

else 

if(gain[3]<0) gain[3]=0; 
if(gain[3]>2) gain[3]=2; 
break; 
case ’n’: 

if((gain[4]>=0) kk (gain[4]<=2)) 

{ 

knobcolor =3* (gainupdown+ 1 ) ; 

_setcolor(knobcolor); 

gainknob=40+(int)(gain[4]*80.0); 

_rectangle(_GFILLINTERIOR,gainknob,448,gainknob+3,468); 

gain[4]=gain[4]+gainupdown*gainstep; 

_ set text posi tion(29,27); 
printf(” %4.3f” ,gain(4] ) ; 

} 

else 

if(gain[4]<0) gain[4]=0; 
if(gain[4]>2) gain[4]=2; 
break; 


} 

} 
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Signal Plotting 


^include <graph.h> 

^include <stdio.h> 

^include <conio.h> 

^include <math.h> 

^include <float.h> 

^include <dos.h> 

^include <bios.h> 

^include <signal.h> 

^include <stdlib.h> 

^include <malloc.h> 

^include ” suspend. h” 

/* Plotdat: plot of signals */ 
void plotdat(void) 

{ 

int plotcenter=148; 
if (plotmode==V) 
kv=plotcenter-(int)(v*10*zoom); 
else 

kv=4+(j+l)*48-(int)(v*2*zoom); 
if(kv>=290) kv=290; 

_setcolor(plotcolor); 

_setpixel(ix,olddat[ix][j4-l]); 

_ setcolor(datcolor ) ; 
_setpixel(ix,kv); 
olddat [ix] [j + 1] = k v ; 

} 



Plots Position Knobs 


^include <graph.h> 
^include <stdio.h> 
^include <conio.h> 
^include <math.h> 
^include <float.h> 
^include <dos.h> 
^include <bios.li> 
^include <signal.h> 
^include <stdlib.h> 
^include <malloc.h> 
^include ”suspend.h” 


/ * positionplot: control and display of position knobs */ 
void positionplot(void) 

{ 

int knobcolor,positionknob,bkcolor; 
knobcolor=10; 
bkcolor=4; 
if(trig==l) 
ix=25; 

sampler=0; /* counter for collecting data in ISR */ 

switcb(ch2) 

{ 


case ’x’: 

_setcolor(0); 

_rectangle(_GFILLINTERIOR,310,320,461,340); 
_ setcolor(bkcolor ) ; 

_rectangle(_GBORDER,310,320,461,340); 

_ setcolor (knobcolor ) ; 

position[0]=position[0]+positionsign*positionstep; 
imposition [0] <- 1 ) position [0] =- 1 ; 



if(position[0]>l) position[0]=l; 
positionknob=385-f (int)(position[0]*75.0); 

_rectangle(_GFILLINTERIOR,positionknob,320,positionknob+l,340); 

break; 

case ’y’ : 

_setcolor(0); 

_rectangle(_GFILLINTERIOR,310,352,461,372); 

_setcolor(bkcolor) ; 

_rectangle(_GBORDER,310,352,461,372); 

_setcolor(knobcolor); 

position [ 1 ] =position [1 ] +positionsign*positionstep; 
if(position[l]<-l) position[l]=-l; 
imposition [1]>1) position[l]=l; 
positionknob=385+(int)(position[l]*75.0); 

_rectangle(_GFILLINTERIOR,positionknob, 352, positionknob+ 1,372); 
break; 

case ’z’: 

_setcolor(0); 

_rectangle(_GFILLINTERIOR,310,384,461,404); 

_ set color ( bkcolor ) ; 

_rectangle(_GBORDER,310,384,461,404); 

_ set color ( knobcolor ) ; 

position[2] =position[2]+positionsign*positionstep; 
if(position[2] <- 1 ) position[2]=- 1 ; 
if(position[2]>l) position[2]=l; 
positionknob=385+(int)(position[2]*75.0); 

_rectangle(_GFILLINTERIOR,positionknob,384,positionknob+l,404); 
break; 
case ’m’: 

_setcolor(0); 

_rectangle(_GFILLINTERIOR,310,416,461,436); 

_ set color ( bkcolor ) ; 

_rectangle(_GBORDER,310,416,461,436); 

_ set color (knobcolor ) ; 

position[3]=position[3]+positionsign*positionstep; 
if(position[3] <- 1 ) position [3] =- 1 ; 



if(position[3]>l) position[3]=l; 

positionknob=385+(int)(position[3]*75.0); 

_rectangle( _ GFILLINTERIOR,positionknob ,416 ,positionknob+ 1 ,436) ; 

break; 
case ’n’: 

_setcolor(0); 

_rectangle(_GFILLINTERIOR,310,448,461,468); 

_setcolor(bkcolor); 

_rectangle(_GBORDER,310,448,461,468); 

_ setcolor (knobcolor) ; 

position [4] =position [4] +positionsign*positionstep ; 
if(position[4] <-l ) position[4]=- 1 ; 
if (positional > 1 ) position[4] = 1 ; 
posi tionknob = 385 + (int ) ( position [4] * 75 . 0 ) ; 

_rectangle(_GFILLINTERIOR,positionknob,448,positionknob+l,468); 

break; 


} 

} 
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Saves Step Responses 


^include <graph.h> 
^include <stdio.h> 
^include <conio.h> 
^include <math.h> 
# include < float. h> 
^include <dos.h> 
^include <bios.h> 
^include <signal.h> 
^include <stdlib.h> 
^include <malloc.h> 
#include ” suspend. h” 


/******** ***************** *******************************************/ 


/* savedat() saves data to file */ 
void savedat(void) 

{ 

char fname[12]; 
int i,j; 

FILE *fp; 

_settextposition(25,64); 
printf(”Save data? y/n\n”); 
ch=(chax)getch(); 


if«ch==’Y’) || (ch==y» 

{ 

_settextposition(25,64); 

printf(” Enter File Name”); 

_settextposition(26,64); 

scanf(”%s”,&fname); 

if((fp=fopen(fname,”w+”))!=NULL) 

{ 


fprintf (fp stepdat = [\n” ) ; 
for(i=l;i<=200;i+-f ) 


{ 
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for(j=0;j<=4;j++) 
fprintf(fp,”%le n ,savdat[j][i]); 
fprintf(fp,”\n”); 

} 

fprintf(fp,”];\n”); 
fprintf(fp,” gains= [\n” ); 

for(i=0;i<=4;i++) 

fprintf(fp,”%le \n”,gain[i]); /* Save the gains */ 

fprintf (fp ] ;\n” ) ; 

fclose(fp); 

} 

else 

printf(”Oops file error”); 

} 


} 



APPENDIX B 


Eddy Current Simulation 


clear 

format short e 

% plant model starting from the analytical model (ls_modl) 
load ls_modl; load varibles; 

% Eddy current effects 
num=[2.1749e-3 1]; 
den=[2.934e-3 1]; 

[ai,bi,ci,di]=tf2ss(num,den); 

[aii ,bii ,cii ,dii] =append( ai ,bi ,ci ,di ,ai ,bi ,ci ,di) ; 

[aiii,biii,ciii,diii]=append(aii,bii,cii,dii,aii,bii,cii,dii); 

[ae,be,ce,de]=append(aiii,biii,ciii,diii,ai,bi,ci,di); 

% adding the dynamics to the plant 
[ap,bp,cp,dp]=series(ae,be,ce,de,A,B,C,D); 


% loading the A B C D matrices of the controller 
load lqg.mat 
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% begin putting the model together 
[ad,bd,cd,dd]=c2dm(ap,bp,cp,dp,T,’zoh’); 

[atotal,btotal,ctotal,dtotal]=series(ac,bc,cc,dc,ad,bd,cd,dd); 

[as,bs,cs,ds]=cloop(atotal,btotal,ctotal,dtotal,-l); 

% Correcting inputs and output of system for position 

bs=bs*p2s; 

cs=s2p*cs; 

ds=s2p*ds*p2s; 

% Performing the simulated step response 

[x,y]=dstep(as,bs,cs,ds,l,40); 

plot(x(:,l)) 
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APPENDIX C 


Numeric Values For Eddy Current Loops 

This list represents the initial startup values calculated from theory. 

Main Coil 

Kz = 8.7736 x 10 5 
Ring 

Lm = 3.873 x 10' 7 Le = 3.8573 x 10' 6 Re = 2.243 x 10' 4 Ke = 4.29 x 10' 6 

Baseplate 

Lm = 4.9056 x 10 ' 5 Le = 5.2523 x 10' 9 Re = 1.4 x 10' 5 Ke=8.508 x 10' 10 

Core 

Lm = 4.9056 x 10' 5 Le = 4.4857 X 10‘ 7 Re = 6 x 10' 3 Ke = 4.601 x 10' 10 
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Appendix D 


OPERA Field Calculating Files 

The data files are generated using the file gradient. comi, which is listed below, in 
OPERA. After the data files are generated the MATLAB script file gradient. m 
is used to find the field gradients. 

GRADIENT.COMI 

grid -#l/2 0 -#1/10 #1/10 #1/10 #1/10 11 1 3 fieldx.dat no 
grid 0 -#l/2 0 #1/10 #1/10 #1/10 1 11 1 fieldy.dat no 
grid 0 0 -#l/2 #1/10 #1/10 #1/10 1 1 11 fieldz.dat no 

GRADIENT.M 

clear 

format long e 
load fieldx.dat 
load fieldy.dat 
load fieldz.dat 
fieldx=fieldx( : ,1 :6) ; 
fieldy=fieldy(:,l:6); 
fieldz=fieldz( : ,1 :6) ; 
bx=fieldy(6,4); 
by=fieldy(6,5); 
bz=fieldy(6,6); 

fieldxl=fieldx([l 4 7 10 13 16 19 22 25 28 31],:); 

fieldx2=fieldx([2 5 8 11 14 17 20 23 26 29 32],:); 

fieldx3=fieldx([3 6 9 12 15 18 21 24 27 30 33],:); 

fieldxl(:,l)=fieldxl(:,l)-fieldxl(6,l); 

fieldx2(:,l)=fieldx2(:,l)-fieldx2(6,l); 

fieldx3(:,l)=fieldx3(:,l)-fieldx3(6,l); 

fieldy ( : ,2 ) =fieldy ( : ,2 )-fieldy (6 ,2) ; 
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fieldz(:,3)=fieldz(:,3)-fieldz(6,3); 

ansi =polyfit (fieldx2( : ,1 ) ,fieldx2( : ,4) ,2) ; 

bxdx=ansl(l,2); 

bxdxdx=2*ansl(l,l); 

ans2 =polyfi t (fieldx2 ( : , 1 ) ,fieldx2 ( : ,5 ) ,2 ) ; 

bydx=ans2(l,2); 

bxdxdy=2*ans2(l,l); 

ans3=polyfit (fieldx2 ( : , 1 ) ,fieldx2( : ,6) ,2) ; 

bxdz=ans3(l,2); 

bxdxdz=2*ans3(l,l); 

ans4=polyfit(fieldy(:,2),fieldy(:,4),2); 

bxdy=ans4(l,2); 

bxdydy=2*ans4(l,l); 

ans5=polyfit(fieldy(:,2),fieldy(:,5),2); 

bydy=ans5(l,2); 

bydydy=2*ans5(l,l); 

ans6=polyfit(fieldy(:,2),fieldy(:,6),2); 

bydz=ans6(l,2); 

bydydz=2*ans6(l,l); 

ans7=polyfit (fieldz( : ,3) ,fieldz(:,4) ,2); 

bzdx=ans7(l,2); 

bxdzdz=2* ans7 (1,1); 

ans8=polyfit (fieldz( : ,3 ) ,fieldz( : ,5) ,2) ; 

bzdy=ans8(l,2); 

bydzdz=2*ans8(l,l); 

ans9=polyfit (fieldz( : ,3 ) ,fieldz( : ,6 ) ,2) ; 

bzdz=ans9(l,2); 

bzdzdz=2*ans9(l,l); 

ans 1 O=polyfit (fieldx 1 ( : , 1 ) ,fieldxl ( : ,5) ,2) ; 

bydxl=ansl0(l,2); 

ansi 1 =polyfit (fieldx3( : , 1 ) ,fieldx3( : ,5) ,2) ; 
bydx3=ansll(l,2); 
l=fieldx2(l,3)-fieldxl(l,3); 
g=[-l , bydxl 

0 , bydx 

1 , bydx3 ]; 



ansl2=polyfit(g(:,l)>6( : ’2),2); 

bxdydz=ansl2(l,2); , 

save data.ma bx by bz bxdx bxdy bxdz bydy bydz bzdz bxdxdx bxdxdy bxdxdz 
bxdydy bxdydz bxdzdz bydydy bydydz bydzdz bzdzdz -ascn 
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